setwd(getwd())
#setwd("D:\thasegawa\ExtremeEvent\R")
library(reshape2)
library(ggplot2)
library(plyr)
library(grid)
library(scales)
library(gdxrrw)
library(magrittr)


regionagmiplabel <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="CAN"] <- "Canada"
  levels(y$Region)[levels(y$Region)=="USA"] <- "USA"
  levels(y$Region)[levels(y$Region)=="BRA"] <- "Brazil"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="IND"] <- "India"
  levels(y$Region)[levels(y$Region)=="OSA"] <- "Other S&C America"
  levels(y$Region)[levels(y$Region)=="FSU"] <- "Former Soviet Union"
  levels(y$Region)[levels(y$Region)=="EUR"] <- "Europe"
  levels(y$Region)[levels(y$Region)=="MEN"] <- "Middle-East/North Afr."
  levels(y$Region)[levels(y$Region)=="SSA"] <- "Sub-Saharan Afr."
  levels(y$Region)[levels(y$Region)=="SEA"] <- "Southeast Asia"
  levels(y$Region)[levels(y$Region)=="OAS"] <- "Other Asia"
  levels(y$Region)[levels(y$Region)=="ANZ"] <- "Australia/New Zealand"
  levels(y$Region)[levels(y$Region)=="NAM"] <- "North America"
  levels(y$Region)[levels(y$Region)=="OAM"] <- "S&C America"
  levels(y$Region)[levels(y$Region)=="AME"] <- "Africa & Middle East"
  levels(y$Region)[levels(y$Region)=="SAS"] <- "Southern Asia"
  z <- y
  return(z)
}
regionagmiplabelshort <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="CAN"] <- "Canada"
  levels(y$Region)[levels(y$Region)=="USA"] <- "USA"
  levels(y$Region)[levels(y$Region)=="BRA"] <- "Brazil"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="IND"] <- "India"
  levels(y$Region)[levels(y$Region)=="OSA"] <- "Other S&C America"
  levels(y$Region)[levels(y$Region)=="FSU"] <- "Former Soviet\nUnion"
  levels(y$Region)[levels(y$Region)=="EUR"] <- "Europe"
  levels(y$Region)[levels(y$Region)=="MEN"] <- "Middle-East/\nNorth Africa"
  levels(y$Region)[levels(y$Region)=="SSA"] <- "Sub-Saharan\nAfrica"
  levels(y$Region)[levels(y$Region)=="SEA"] <- "Southeast\nAsia"
  levels(y$Region)[levels(y$Region)=="OAS"] <- "Other\nAsia"
  levels(y$Region)[levels(y$Region)=="ANZ"] <- "Australia/\nNew Zealand"
  levels(y$Region)[levels(y$Region)=="NAM"] <- "North\nAmerica"
  levels(y$Region)[levels(y$Region)=="OAM"] <- "S&C\n America"
  levels(y$Region)[levels(y$Region)=="AME"] <- "Africa & Middle East"
  levels(y$Region)[levels(y$Region)=="SAS"] <- "Southern Asia"
  z <- y
  return(z)
}


fregionagmipnameshort <- function(xx){
  if (xx=="WLD"){
    outname <- paste("world")
  } else if (xx=="CAN"){
    outname <- paste("Canada")
  } else if (xx=="USA"){
    outname <- paste("USA")
  } else if (xx=="BRA"){
    outname <- paste("Brazil")
  } else if (xx=="CHN"){
    outname <- paste("China")
  } else if (xx=="IND"){
    outname <- paste("India")
  } else if (xx=="OSA"){
    outname <- paste("Other S&C\nAmerica")
  } else if (xx=="FSU"){
    outname <- paste("Former Soviet\nUnion")
  } else if (xx=="EUR"){
    outname <- paste("Europe")
  } else if (xx=="MEN"){
    outname <- paste("Middle-East/North Afr.")
  } else if (xx=="SSA"){
    outname <- paste("Sub-Saharan\nAfrica")
  } else if (xx=="SEA"){
    outname <- paste("Southeast\nAsia")
  } else if (xx=="OAS"){
    outname <- paste("Other Asia")
  } else if (xx=="ANZ"){
    outname <- paste("Australia/New Zealand")
  } else if (xx=="NAM"){
    outname <- paste("North\nAmerica")
  } else if (xx=="OAM"){
    outname <- paste("S&C\nAmerica")
  } else if (xx=="AME"){
    outname <- paste("Africa &\nMiddle East")
  } else if (xx=="SAS"){
    outname <- paste("Southern\nAsia")
  } else {
    outname <- c("aaaaaaaa")
  }
  return(outname)
}

fregionagmipname <- function(xx){
  if (xx=="WLD"){
    outname <- paste("world")
  } else if (xx=="CAN"){
    outname <- paste("Canada")
  } else if (xx=="USA"){
    outname <- paste("USA")
  } else if (xx=="BRA"){
    outname <- paste("Brazil")
  } else if (xx=="CHN"){
    outname <- paste("China")
  } else if (xx=="IND"){
    outname <- paste("India")
  } else if (xx=="OSA"){
    outname <- paste("Other S&C America")
  } else if (xx=="FSU"){
    outname <- paste("Former Soviet Union")
  } else if (xx=="EUR"){
    outname <- paste("Europe")
  } else if (xx=="MEN"){
    outname <- paste("Middle-East/North Afr.")
  } else if (xx=="SSA"){
    outname <- paste("Sub-Saharan Afr.")
  } else if (xx=="SEA"){
    outname <- paste("Southeast Asia")
  } else if (xx=="OAS"){
    outname <- paste("Other Asia")
  } else if (xx=="ANZ"){
    outname <- paste("Australia/New Zealand")
  } else if (xx=="NAM"){
    outname <- paste("North America")
  } else if (xx=="OAM"){
    outname <- paste("S&C America")
  } else if (xx=="AME"){
    outname <- paste("Africa & Middle East")
  } else if (xx=="SAS"){
    outname <- paste("Southern Asia")
  } else {
    outname <- c("aaaaaaaa")
  }
  return(outname)
}


regionlabel <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="USA"] <- "United States"
  levels(y$Region)[levels(y$Region)=="XE25"] <- "EU25"
  levels(y$Region)[levels(y$Region)=="XER"] <- "Rest of Europe"
  levels(y$Region)[levels(y$Region)=="TUR"] <- "Turkey"
  levels(y$Region)[levels(y$Region)=="XOC"] <- "Oceania"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="IND"] <- "India"
  levels(y$Region)[levels(y$Region)=="JPN"] <- "Japan"
  levels(y$Region)[levels(y$Region)=="XSE"] <- "Southeast Asia"
  levels(y$Region)[levels(y$Region)=="XSA"] <- "Rest of South Asia"
  levels(y$Region)[levels(y$Region)=="CAN"] <- "Canada"
  levels(y$Region)[levels(y$Region)=="BRA"] <- "Brazil"
  levels(y$Region)[levels(y$Region)=="XLM"] <- "Rest of South America"
  levels(y$Region)[levels(y$Region)=="CIS"] <- "CIS"
  levels(y$Region)[levels(y$Region)=="XME"] <- "Middle East"
  levels(y$Region)[levels(y$Region)=="XNF"] <- "North Africa"
  levels(y$Region)[levels(y$Region)=="XAF"] <- "Sub-Sa Africa"
z <- y
return(z)
}

regionlabelagg <- function(y){
  levels(y$Region)[levels(y$Region)=="WLD"] <- "World"
  levels(y$Region)[levels(y$Region)=="US"] <- "United States"
  levels(y$Region)[levels(y$Region)=="LAM"] <- "Latin America"
  levels(y$Region)[levels(y$Region)=="EU"] <- "EU"
  levels(y$Region)[levels(y$Region)=="CHN"] <- "China"
  levels(y$Region)[levels(y$Region)=="SEAsia"] <- "Southeast Asia"
  levels(y$Region)[levels(y$Region)=="SAsia"] <- "South Asia"
  levels(y$Region)[levels(y$Region)=="Afr_ME"] <- "Africa & Middle East"
  levels(y$Region)[levels(y$Region)=="FSU"] <- "CIS"
  levels(y$Region)[levels(y$Region)=="Other"] <- "Other"
  z <- y
  return(z)
}
croplabel <- function(y){
  levels(y$Crop)[levels(y$Crop)=="PDR"] <- "rice"
  levels(y$Crop)[levels(y$Crop)=="WHT"] <- "wheat"
  levels(y$Crop)[levels(y$Crop)=="GRO"] <- "other grains"
  levels(y$Crop)[levels(y$Crop)=="OSD"] <- "oil crops"
  z <- y
  return(z)
}

varlabel <- function(y){
  levels(y$VAR)[levels(y$VAR)=="POPU"] <- "Population at risk of hunger \n [million]"
  levels(y$VAR)[levels(y$VAR)=="CALO"] <- "Per capita food consumption \n [kcal/person/day]"
  levels(y$VAR)[levels(y$VAR)=="YEXO"] <- "Yield \n [2005=1]"
  levels(y$VAR)[levels(y$VAR)=="XPRP"] <- "Price index \n [base year=1]"
  levels(y$VAR)[levels(y$VAR)=="AREA"] <- "Cropland area \n [Mha]"
  levels(y$VAR)[levels(y$VAR)=="PROD"] <- "Production \n [million tonne]"
  
  z <- y
  return(z)
}

varlabelnounit <- function(y){
  levels(y$VAR)[levels(y$VAR)=="POPU"] <- "Population at risk of hunger"
  levels(y$VAR)[levels(y$VAR)=="CALO"] <- "Per capita food consumption"
  levels(y$VAR)[levels(y$VAR)=="YEXO"] <- "Yield"
  levels(y$VAR)[levels(y$VAR)=="XPRP"] <- "Price index"
  levels(y$VAR)[levels(y$VAR)=="AREA"] <- "Cropland area"
  levels(y$VAR)[levels(y$VAR)=="PROD"] <- "Production"
  
  z <- y
  return(z)
}



fsize <- 13
MyThemeLine <- theme_grey() +
  theme(
    panel.border=element_rect(fill=NA),
    title=element_text(size=fsize,face="bold"),
    axis.title=element_text(size=fsize,face="bold"),
    axis.text = element_text(hjust=1,size = fsize, angle = 0),
    axis.line=element_line(colour="black"),
    panel.background=element_rect(fill = "white"),
    #    panel.background=element_blank(),
    #    panel.grid.major=element_line(linetype="dashed",colour="grey",size=0.5),
    panel.grid.major=element_blank(),
    strip.background=element_rect(fill="white", colour="white"),
    strip.text.x = element_text(size=fsize, colour = "black", angle = 0,face="bold"),
    legend.title = element_text(size=fsize),
    legend.text = element_text(size=fsize)
  )

#Defining function
#Making two figures
fourfigure <- function(p1,p2,p3,p4,x){
  grid.newpage()
  l <- grid.layout( nrow=2, ncol=2, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1, 1),c("null","null")))
  png(filename=x, width=2400, height=1500,res=200)
  v <- viewport(layout = l) # [l]を作図領域として[v]に格納す???
  pushViewport(v) #グラフを埋め込む前に???pushViewport()関数に作図領域を指定する???
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p3, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  popViewport()  #作図を終???する
  dev.off()  
  return()
}
nc <- 0
xhist <- 0
Dataextract <- function(x,y,v,reg,cr,yr){
  xhist <- x[x$VAR==v & x$Year==yr & x$Region==reg & x$SSP=="SSP2" & x$Crop==cr, c(-1,-2,-3,-4,-10)]
  nocc1 <- y[y$VAR==v & y$Year==yr & y$Region==reg & y$SSP=="SSP2" & y$Crop==cr, c(6)]  
  names(xhist) <- c("GCM","RCP","CO2fert","SYR","N","Value")	#Rename the column
  minv <- min(min(xhist[6]),nocc1)
  if (varname=="PPOP_UNSH"){
  maxv <- max(max(xhist[6]),nocc1)
    bw <- maxv/100
  } else {
    maxv <- max(max(xhist[6]),nocc1)
    bw <- diff(range(xhist$Value))/100
  }
  
  mainlabel <- paste(reg,cr,yr)
  interact<-interaction(xhist$RCP, xhist$CO2fert, sep=" : ")
  
  mu <- ddply(xhist, "interact", summarise, grp.mean=mean(Value))
  v5 <- ddply(xhist,"interact",summarise, grp.quantile=quantile(Value,0.05))
  plotx <- ggplot() + 
    geom_density(data=xhist, aes(x=Value,fill=interact,y=..density..,color=interact),alpha = 0,size=2, position = "identity") +  
    geom_histogram(data=xhist, aes(x=Value,fill=interact,y=..density..,color=interact),alpha = 0.3, position = "identity") +
    #    geom_histogram(alpha = 0.3, position = "identity",  binwidth=bw, color="black") +
    xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
    ggtitle(mainlabel) +
    ylab("density") + xlab(labelname) +
    MyThemeLine +
    geom_vline(data=mu,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1) +
  #  geom_vline(data=v5,aes(xintercept=grp.quantile, color=interact), linetype="solid", size=1) +
    #  geom_text(data=mu,aes(x=grp.mean,y=0.023,label=grp.mean),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
    # geom_text(data=v5,aes(x=grp.quantile,y=0.023,label=grp.quantile),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001) +
  geom_vline(aes(xintercept=nocc1), color="black", size=1.2) 
  #  geom_text(aes(x=nocc1,y=0.023,label=nocc1),angle=-90,size=5,vjust=-bw/10,hjust=-0.0001)
  #  geom_area(mapping = aex(y = ifelse(y<0.2)))
#  ggsave(plot1,filename=filenamehist, dpi=250,width=8, height=6)
  return(plotx)  
}

twofigure <- function(p1,p2,x){
  grid.newpage()
  l <- grid.layout( nrow=1, ncol=2, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1, 1),c("null","null")))
  png(filename=x, width=2400, height=1500,res=200)
  v <- viewport(layout = l) # [l]を作図領域として[v]に格納す???
  pushViewport(v) #グラフを埋め込む前に???pushViewport()関数に作図領域を指定する???
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  popViewport()  #作図を終???する
  dev.off()  
  return()
}

#Making five figures
fivefigure <- function(p1,p2,p3,p4,p5,x,l,wd,ht,rs,x2){
  grid.newpage()
  png(filename=x, width=wd, height=ht,res=rs)
  v <- viewport(layout = l)
  pushViewport(v)
  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=3,layout.pos.col=1))
  print(p4, vp = viewport(layout.pos.row=3,layout.pos.col=2))
  print(p5, vp = viewport(layout.pos.row=4,layout.pos.col=1))
  #  print(title=x2)
  popViewport()
  dev.off()  
  return()
}


#---- making directory
dirname0 <- paste('../output/png/', sep = '')
dir.create(dirname0)

var.in <- c("PCAL","PPOP_UNSH","YEXO","XPRP")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  #  dirname <- paste('../output/png/', varname, sep = '')
  dirname <- paste(dirname0, varname, sep = '')
  dir.create(dirname)
  
}

linepalette <- c("#4DAF4A","#377EB8","#E41A1C","#FF7F00","#984EA3","#FFFF33","#A65628","#F781BF","#8DD3C7","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F")

#--- data load

# extreme diff from nocc
colnames <- c("Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","VAR","Value")
pdata_diff <- rgdx.param('../data/Pagmipout.gdx','Pagmipout_diff')
names(pdata_diff) <- colnames

# extreme absolute
pdata <- rgdx.param('../data/Pagmipout.gdx','Pagmipout')
names(pdata) <- colnames
pdata <- pdata[pdata$N!="N0", c(1,2,3,4,5,6,7,8,9,10,11)]

# mean
colnames_long <- c("Year","Region","SSP","GCM","RCP","CO2fert","Crop","VAR","Value")
pdata_long <- rgdx.param('../data/Pagmipout_rev.gdx','Pagmipout_longterm')
names(pdata_long) <- colnames_long

pdata$RCP <- factor(pdata$RCP, levels=c("RCP85","RCP26"))
pdata$Region <- factor(pdata$Region, levels=c("WLD","CAN","USA","BRA","CHN","IND","OSA","FSU","EUR","MEN","SSA","SEA","OAS","ANZ","NAM","OAM","AME","SAS"))

#--- data load end

col<- c('red','blue','cyan','orange')

#---

#1 Boxplot
y <- pdata_diff[pdata_diff$SSP=="SSP2" & pdata_diff$Year==2050 & pdata_diff$CO2fert=="noco2" & pdata_diff$Region=="WLD", c(-1,-2,-6)]
y <- y[(y$Crop=="AGR"& y$VAR=="POPU") | (y$Crop=="AGR"& y$VAR=="CALO") | (y$Crop=="CR5" & y$VAR=="YEXO") | (y$Crop=="AGR" & y$VAR=="XPRP") | (y$Crop=="CRP" & y$VAR=="AREA") | (y$Crop=="CR5" & y$VAR=="PROD") ,c(1,2,3,4,5,6,7,8)]

interact<-interaction(y$RCP, y$GCM, sep=" : ")
y <- varlabel(y)

gbox <- ggplot() + 
  geom_boxplot(data=y, aes(x=RCP, y=Value, group=RCP),fill="blue", alpha=0.3, colour="Blue") +
  geom_boxplot(data=y, aes(x=interact, y=Value, group=interact),fill="blue", alpha=0.3, colour="Blue") +
#  geom_jitter(data=y, aes(x=interact, y=Value, group=RCP)) +
  ylab("") + xlab("") +	
  MyThemeLine +
#  geom_point(data=z, aes(x=Region, y=NoCC),size=10,shape=73,colour="black")
  theme(
    axis.text = element_text(hjust=1,size = fsize, angle = 90)
  )

gbox2 <- gbox + facet_wrap(~VAR,ncol=3,scales = 'free_y') 
plot(gbox2)

filenamebox <- paste('../output/png/boxplot_gcm_2050_N100rev.png', sep = '')
ggsave(gbox2,filename=filenamebox, dpi=250,width=12, height=8)



x <- pdata_diff[pdata_diff$SSP=="SSP2" & pdata_diff$CO2fert=="noco2" & pdata_diff$Region=="WLD", c(-2,-6)]
x <- x[(x$Crop=="AGR"& x$VAR=="POPU") | (x$Crop=="AGR"& x$VAR=="CALO") | (x$Crop=="CR5" & x$VAR=="YEXO") | (x$Crop=="AGR"& x$VAR=="XPRP") | (x$Crop=="CRP"& x$VAR=="AREA") | (x$Crop=="CR5"& x$VAR=="PROD"),c(1,2,3,4,5,6,7,8,9)]
x <- varlabel(x)

interact2<-interaction(x$RCP, x$Year, sep=" : ")

gbox3 <- ggplot() + 
#  geom_jitter(data=x, aes(x=interact2, y=Value), size = 1, colour='gray50') +
  geom_boxplot(data=x, aes(x=interact2, y=Value, group = interact2), fill="blue", alpha=0.3, colour="Blue",outlier.shape = NA) +
  ylab("") + xlab("") +	
  MyThemeLine +
  #  geom_point(data=z, aes(x=Region, y=NoCC),size=10,shape=73,colour="black")
  theme(
    axis.text = element_text(hjust=1, size = fsize, angle = 90)
  )

gbox4 <- gbox3 + facet_wrap(~VAR,ncol=3,scales = 'free_y') 
plot(gbox4)

filenamebox2 <- paste('../output/png/boxplot_year_N100rev.png', sep = '')
ggsave(gbox4,filename=filenamebox2, dpi=250,width=10, height=8)

# F test 

z <- pdata[pdata$SSP=="SSP2" & (pdata$Year==2050 | pdata$Year==2020) & pdata$CO2fert=="noco2" & pdata$Region=="WLD", c(-2,-6)]
z0 <- z[z$Crop=="AGR" & z$VAR=="POPU",c(-7,-8)]

zvar0 <- ddply(z0,c("RCP","Year"), summarise, variance = var(Value))

z1 <- ddply(z0,c("RCP","Year"), transform, mean = mean(Value))
z1[7] <- z1[7]/z1[8]

zvar <- ddply(z1,c("RCP","Year"), summarise, variance = var(Value))

data2026 <- z1[z1$Year=="2020" & z1$RCP=="RCP26",c(7)] 
data5026 <- z1[z1$Year=="2050" & z1$RCP=="RCP26",c(7)] 
data2085 <- z1[z1$Year=="2020" & z1$RCP=="RCP85",c(7)] 
data5085 <- z1[z1$Year=="2050" & z1$RCP=="RCP85",c(7)] 

library(tidyr)

var.test(data2026, data5026)
var.test(data2085, data5085)

# Boxplot 3 -- decomposition analysis

z <- pdata[pdata$SSP=="SSP2" & (pdata$Year==2050 | pdata$Year==2020) & pdata$CO2fert=="noco2" & pdata$Region=="WLD", c(-2,-6)]
z <- z[(z$Crop=="AGR" & z$VAR=="POPU") | (z$Crop=="AGR"& z$VAR=="CALO") | (z$Crop=="CR5" & z$VAR=="YEXO") ,c(1,2,3,4,5,6,7,8,9)]

interact<-interaction(z$RCP, z$GCM, sep=" : ")
z1 <- varlabel(z)

gboxz <- ggplot() + 
  geom_boxplot(data=z1, aes(x=RCP, y=Value, group=RCP),fill="blue", alpha=0.3, colour="Blue") +
  geom_boxplot(data=z1, aes(x=interact, y=Value, group=interact),fill="blue", alpha=0.3, colour="Blue") +
  #  geom_jitter(data=y, aes(x=interact, y=Value, group=RCP)) +
  ylab("") + xlab("") +	
  MyThemeLine +
  #  geom_point(data=z, aes(x=Region, y=NoCC),size=10,shape=73,colour="black")
  theme(
    axis.text = element_text(hjust=1,size = fsize, angle = 90)
  )

gboxz2 <- gboxz + facet_grid(VAR~Year,scales = 'free_y') 
plot(gboxz2)

filenameboxz <- paste('../output/png/boxplot_gcm_year_rcprev.png', sep = '')
ggsave(gboxz2,filename=filenameboxz, dpi=250,width=7, height=9)


zmean <- ddply(z,c("VAR","RCP","GCM","Year"), summarise, mean = mean(Value))
zsd <- ddply(z,c("VAR","RCP","GCM","Year"), summarise, sd = sd(Value))
zmean2 <- ddply(z,c("VAR","RCP","Year"), summarise, mean = mean(Value))
zsd2 <- ddply(z,c("VAR","RCP","Year"), summarise, sd = sd(Value))

zstat_gcm <- cbind(zmean,zsd)
zstat_rcp <- cbind(zmean2,zsd2)
zstat_gcm[11] <- zstat_gcm[10]/zstat_gcm[5]
zstat_rcp[9] <- zstat_rcp[8]/zstat_rcp[4]

zstat_gcm2 <- zstat_gcm[c(1,2,3,4,5,10,11)]
zstat_rcp2 <- zstat_rcp[c(1,2,3,4,8,9)]

colnamesz_gcm <- c("VAR","RCP","GCM","Year","Mean","SD","CV")
colnamesz_rcp <- c("VAR","RCP","Year","Mean","SD","CV")

names(zstat_gcm2) <- colnamesz_gcm
names(zstat_rcp2) <- colnamesz_rcp

zstat_gcm3 <- varlabelnounit(zstat_gcm2)

gbarz <- ggplot() + 
#  geom_bar(data=zstat_gcm3, aes(x=GCM, y=CV,fill=Year), position="dodge", stat="identity") +
  geom_boxplot(data=zstat_gcm3, aes(x=Year, y=CV, group=Year),fill="blue", alpha=0.3, colour="Blue") +
  ylab("CV") + xlab("") +	
  MyThemeLine +
  theme(
    axis.text = element_text(hjust=0.5)
  )

gbarz2 <- gbarz + facet_grid(VAR~RCP,scales = 'free_y') 
plot(gbarz2)

#filenamebarz <- paste('../output/png/bar_gcm_year_rcp_cv.png', sep = '')
filenamebarz <- paste('../output/png/box_gcm_year_rcp_cv.png', sep = '')
ggsave(gbarz2,filename=filenamebarz, dpi=250,width=5, height=7)

zstat_gcm4 <- varlabel(zstat_gcm2)

gbarrcpz <- ggplot() + 
#  geom_bar(data=zstat_gcm4, aes(x=GCM, y=Mean,fill=Year), position="dodge", stat="identity") +
  geom_boxplot(data=zstat_gcm4, aes(x=Year, y=Mean, group=Year),fill="blue", alpha=0.3, colour="Blue") +
  ylab("Mean") + xlab("") +	
  MyThemeLine +
  theme(
    axis.text = element_text(hjust=0.5,size = fsize)
  )

gbarrcpz2 <- gbarrcpz + facet_grid(VAR~RCP, scales = 'free_y') 
plot(gbarrcpz2)

#filenamebarzmean <- paste('../output/png/bar_gcm_year_rcp_mean.png', sep = '')
filenamebarzmean <- paste('../output/png/box_gcm_year_rcp_mean.png', sep = '')
ggsave(gbarrcpz2,filename=filenamebarzmean, dpi=250,width=5, height=7)


#----- #1-2 Boxplot across regions ----

#s <- pdata[pdata$Year==yearhist & pdata$Crop=="AGR" & pdata$CO2fert=="noco2", c(-1,-6,-9)]
#s <- s[(s$VAR=="POPU" | s$VAR=="CALO"),c(1,2,3,4,5,6,7,8)]

s <- pdata[pdata$SSP=="SSP2" & pdata$Year==yearhist & pdata$CO2fert=="noco2", c(-1,-6)]
s <- s[(s$Crop=="AGR"& s$VAR=="POPU") | (s$Crop=="AGR"& s$VAR=="CALO") | (s$Crop=="CR5" & s$VAR=="YEXO") | (s$Crop=="AGR"& s$VAR=="XPRP"), c(1,2,3,4,5,6,7,8,9)]
s <- s[s$Region=="CHN" | s$Region=="IND" | s$Region=="FSU" | s$Region=="MEN" | 
       s$Region=="SSA" | s$Region=="SEA" | s$Region=="OAS" | s$Region=="OAM" |
       s$Region=="EUR" | s$Region=="ANZ" | s$Region=="NAM" 
       ,c(1,2,3,4,5,6,7,8,9)]
s$VAR <- factor(s$VAR, levels=c("YEXO","XPRP","CALO","POPU"))
s$Region <- factor(s$Region, levels=c("SSA","IND","CHN","SEA","OAS","FSU","MEN","OAM","EUR","ANZ","NAM"))
s <- regionagmiplabelshort(s)
s <- varlabel(s)

mins <- ddply(s,c("Region","RCP","VAR"),summarize,min=min(Value))
maxs <- ddply(s,c("Region","RCP","VAR"),summarize,max=max(Value))
v99 <- ddply(s,c("Region","RCP","VAR"),summarize,v99=quantile(Value,0.99))
v01 <- ddply(s,c("Region","RCP","VAR"),summarize,v01=quantile(Value,0.01)) 
stats <- cbind(regmins,regmaxs,regv01,regv99)

nocc2 <- pdata_long[pdata_long$Year==yearhist & pdata_long$SSP=="SSP2" & pdata_long$GCM=="nocc", c(-1,-3,-4,-5,-6)]
#nocc2 <- nocc2[(nocc2$VAR=="POPU" | nocc2$VAR=="CALO" | nocc2$VAR=="XPRP"),c(1,2,3)]
nocc2 <- nocc2[(nocc2$Crop=="AGR"& nocc2$VAR=="POPU") | (nocc2$Crop=="AGR"& nocc2$VAR=="CALO") | (nocc2$Crop=="CR5" & nocc2$VAR=="YEXO") | (nocc2$Crop=="AGR"& nocc2$VAR=="XPRP"),c(-2)]
nocc2 <- nocc2[nocc2$Region=="CHN" | nocc2$Region=="IND" | nocc2$Region=="FSU" | nocc2$Region=="MEN" | 
                 nocc2$Region=="SSA" | nocc2$Region=="SEA" | nocc2$Region=="OAS" | nocc2$Region=="OAM"|
                 nocc2$Region=="EUR" | nocc2$Region=="ANZ" | nocc2$Region=="NAM" 
               ,c(1,2,3)]
nocc2 <- regionagmiplabelshort(nocc2)
nocc2 <- varlabel(nocc2)

gboxreg <- ggplot() + 
  geom_boxplot(data=s, aes(x=RCP, y=Value, fill=RCP)) +
  geom_boxplot(data=s, aes(x=RCP, y=Value, fill=RCP)) +
  scale_fill_manual(name="RCP",values=col) +
  ylab("") +
  xlab("") +	
  MyThemeLine +
  theme(
    axis.text.x = element_text(hjust=1,size = fsize, angle = 90),
    legend.position = "none",
    strip.text = element_text(size=fsize, colour = "black", angle = 0,face="plain")
  ) +
 geom_hline(data=nocc2,aes(yintercept=Value), color="black", size=1)
  
gboxreg2 <- gboxreg + facet_grid(VAR~Region,scales = 'free')
plot(gboxreg2)

filenamebox <- paste('../output/png/boxplot2_2050_region_withhunger_4index_N100.png', sep = '')
filenamebox <- paste('../output/png/boxplot2_2050_allreg_4index_N100.png', sep = '')
ggsave(gboxreg2,filename=filenamebox, dpi=250,width=18, height=12)


#----End --- Boxplot across regions----

#----- #1-2 Boxplot across regions with co2 and noco2 ----
rcphist <- c("RCP85")
col2<- c('red','orange')

rcphist <- c("RCP26")
col2<- c('blue','cyan')

k <- pdata[pdata$SSP=="SSP2" & pdata$Year==yearhist & pdata$RCP==rcphist, c(-1,-5)]
k <- k[(k$Crop=="AGR"& k$VAR=="POPU") | (k$Crop=="AGR"& k$VAR=="CALO") | (k$Crop=="CR5" & k$VAR=="YEXO") | (k$Crop=="AGR"& k$VAR=="XPRP"), c(1,2,3,4,5,6,7,8,9)]
k <- k[k$Region=="CHN" | k$Region=="IND" | k$Region=="FSU" | k$Region=="MEN" | 
       k$Region=="SSA" | k$Region=="SEA" | k$Region=="OAS" | k$Region=="OAM" 
#       k$Region=="EUR" | k$Region=="ANZ" | k$Region=="NAM" 
       ,c(1,2,3,4,5,6,7,8,9)]
k$VAR <- factor(k$VAR, levels=c("YEXO","XPRP","CALO","POPU"))
k$Region <- factor(k$Region, levels=c("SSA","IND","CHN","SEA","OAS","FSU","MEN","OAM","EUR","ANZ","NAM"))
k$CO2fert <- factor(k$CO2fert, levels=c("noco2","co2"))
k <- regionagmiplabelshort(k)
k <- varlabel(k)

mink <- ddply(k,c("Region","CO2fert","VAR"),summarize,min=min(Value))
maxk <- ddply(k,c("Region","CO2fert","VAR"),summarize,max=max(Value))

statk <- cbind(mink,maxk)

nocc2 <- pdata_long[pdata_long$Year==yearhist & pdata_long$SSP=="SSP2" & pdata_long$GCM=="nocc", c(-1,-3,-4,-5,-6)]
#nocc2 <- nocc2[(nocc2$VAR=="POPU" | nocc2$VAR=="CALO" | nocc2$VAR=="XPRP"),c(1,2,3)]
nocc2 <- nocc2[(nocc2$Crop=="AGR"& nocc2$VAR=="POPU") | (nocc2$Crop=="AGR"& nocc2$VAR=="CALO") | (nocc2$Crop=="CR5" & nocc2$VAR=="YEXO") | (nocc2$Crop=="AGR"& nocc2$VAR=="XPRP"),c(-2)]
nocc2 <- nocc2[nocc2$Region=="CHN" | nocc2$Region=="IND" | nocc2$Region=="FSU" | nocc2$Region=="MEN" | 
               nocc2$Region=="SSA" | nocc2$Region=="SEA" | nocc2$Region=="OAS" | nocc2$Region=="OAM"
#               nocc2$Region=="EUR" | nocc2$Region=="ANZ" | nocc2$Region=="NAM" 
               ,c(1,2,3)]
nocc2 <- regionagmiplabelshort(nocc2)
nocc2 <- varlabel(nocc2)

gboxregco2 <- ggplot() + 
  geom_boxplot(data=k, aes(x=CO2fert, y=Value, fill=CO2fert)) +
  geom_boxplot(data=k, aes(x=CO2fert, y=Value, fill=CO2fert)) +
  scale_fill_manual(name="CO2fert",values=col2) +
  ylab("") +
  xlab("") +	
  MyThemeLine +
  theme(
    axis.text.x = element_text(hjust=1,size = fsize, angle = 45),
    legend.position = "none",
    strip.text = element_text(size=fsize, colour = "black", angle = 0,face="plain")
  ) +
  geom_hline(data=nocc2,aes(yintercept=Value), color="black", size=1)

gboxregco2_2 <- gboxregco2 + facet_grid(VAR~Region,scales = 'free')
plot(gboxregco2_2)

#filenameboxco2 <- paste('../output/png/boxplot2_2050_region_withhunger_4index_co2fert_',rcphist,'_N100.png', sep = '')
filenameboxco2 <- paste('../output/png/boxplot2_2050_allreg_4index_co2fert_',rcphist,'_N100.png', sep = '')
ggsave(gboxregco2_2,filename=filenameboxco2, dpi=250,width=12, height=12)


#----End --- Boxplot across regions with co2 and noco2----


#----- #1-2 Boxplot across regions  with co2 and noco2 and RCP2.6 and RCP8.5 ----
col3<- c('red','blue','orange','cyan')

yearhist <- c("2050")

s <- pdata[pdata$SSP=="SSP2" & pdata$Year==yearhist, c(-1)]
s <- s[(s$Crop=="AGR"& s$VAR=="POPU") | (s$Crop=="AGR"& s$VAR=="CALO") | (s$Crop=="CR5" & s$VAR=="YEXO") | (s$Crop=="AGR"& s$VAR=="XPRP"), c(1,2,3,4,5,6,7,8,9,10)]
s <- s[s$Region=="WLD" | s$Region=="CHN" | s$Region=="IND" | s$Region=="FSU" | s$Region=="MEN" | 
         s$Region=="SSA" | s$Region=="SEA" | s$Region=="OAS" | s$Region=="OAM"
#         s$Region=="EUR" | s$Region=="ANZ" | s$Region=="NAM" 
       ,c(1,2,3,4,5,6,7,8,9,10)]
s$VAR <- factor(s$VAR, levels=c("YEXO","XPRP","CALO","POPU"))
s$Region <- factor(s$Region, levels=c("WLD","SSA","IND","CHN","SEA","OAS","FSU","MEN","OAM","EUR","ANZ","NAM"))
s <- regionagmiplabelshort(s)
s <- varlabel(s)

s$RCP <- factor(s$RCP, levels=c("RCP85","RCP26"))
s$CO2fert <- factor(s$CO2fert, levels=c("noco2","co2"))
interact<-interaction(s$RCP, s$CO2fert, sep="_")

mins <- ddply(s,c("Region","interact","VAR"),summarize,min=min(Value))
maxs <- ddply(s,c("Region","interact","VAR"),summarize,max=max(Value))
v99 <- ddply(s,c("Region","interact","VAR"),summarize,v99=quantile(Value,0.99))
v01 <- ddply(s,c("Region","interact","VAR"),summarize,v01=quantile(Value,0.01)) 
v05 <- ddply(s,c("Region","interact","VAR"),summarize,v50=mean(Value))
stats <- cbind(mins,maxs,v01,v99,v05)
stats2 <- stats[c(1,2,3,4,8,12,16,20)]

nocc2 <- pdata_long[pdata_long$Year==yearhist & pdata_long$SSP=="SSP2" & pdata_long$GCM=="nocc", c(-1,-3,-4,-5,-6)]
#nocc2 <- nocc2[(nocc2$VAR=="POPU" | nocc2$VAR=="CALO" | nocc2$VAR=="XPRP"),c(1,2,3)]
nocc2 <- nocc2[(nocc2$Crop=="AGR"& nocc2$VAR=="POPU") | (nocc2$Crop=="AGR"& nocc2$VAR=="CALO") | (nocc2$Crop=="CR5" & nocc2$VAR=="YEXO") | (nocc2$Crop=="AGR"& nocc2$VAR=="XPRP"),c(-2)]
nocc2 <- nocc2[nocc2$Region=="WLD" | nocc2$Region=="CHN" | nocc2$Region=="IND" | nocc2$Region=="FSU" | nocc2$Region=="MEN" | 
                 nocc2$Region=="SSA" | nocc2$Region=="SEA" | nocc2$Region=="OAS" | nocc2$Region=="OAM"
#                 nocc2$Region=="EUR" | nocc2$Region=="ANZ" | nocc2$Region=="NAM" 
               ,c(1,2,3)]
nocc2 <- regionagmiplabelshort(nocc2)
nocc2 <- varlabel(nocc2)

gboxregco2rcp <- ggplot() + 
  geom_boxplot(data=s, aes(x=interact, y=Value, fill=interact)) +
  scale_fill_manual(name="RCP_CO2fert",values=col3) +
  ylab("") +
  xlab("") +	
  MyThemeLine +
  theme(
    axis.text.x = element_text(hjust=1,size = fsize, angle = 90),
    legend.position = "none",
    strip.text = element_text(size=fsize, colour = "black", angle = 0,face="plain")
  ) +
  geom_hline(data=nocc2,aes(yintercept=Value), color="black", size=1)

gboxregco2rcp2<- gboxregco2rcp + facet_grid(VAR~Region,scales = 'free')
plot(gboxregco2rcp2)


filenameboxco2rcp <- paste('../output/png/boxplot2_2050_region_withhunger_4index_N100_co2rcprev.pdf', sep = '')
ggsave(gboxregco2rcp2,filename=filenameboxco2rcp, dpi=250,width=13, height=12)
filenameboxco2rcp <- paste('../output/pdf/Figure3.pdf', sep = '')
ggsave(gboxregco2rcp2,filename=filenameboxco2rcp, dpi=250,width=13, height=12)

filenameboxco2rcp <- paste('../output/png/boxplot2_2050_allreg_4index_N100_co2rcprev.png', sep = '')
ggsave(gboxregco2rcp2,filename=filenameboxco2rcp, dpi=250,width=18, height=12)


#----End --- Boxplot across regions  with co2 and noco2 and RCP2.6 and RCP8.5----

#----- #1-2 Boxplot across regions  with SSP2 and SSP3 and RCP2.6 and RCP8.5 ----
col3.2<- c('red','violet','blue','lightblue')
col3.3<- c('black','grey')

yearhist <- c("2050")

s <- pdata[(pdata$SSP=="SSP2" | pdata$SSP=="SSP3") & pdata$CO2fert=="noco2" & pdata$Year==yearhist, c(-1,-6)]
s <- s[(s$Crop=="AGR"& s$VAR=="POPU") | (s$Crop=="AGR"& s$VAR=="CALO") | (s$Crop=="CR5" & s$VAR=="YEXO") | (s$Crop=="AGR"& s$VAR=="XPRP"), c(1,2,3,4,5,6,7,8,9)]
s <- s[s$Region=="WLD" | s$Region=="CHN" | s$Region=="IND" | s$Region=="FSU" | s$Region=="MEN" | 
         s$Region=="SSA" | s$Region=="SEA" | s$Region=="OAS" | s$Region=="OAM" |
         s$Region=="EUR" | s$Region=="ANZ" | s$Region=="NAM" 
       ,c(1,2,3,4,5,6,7,8,9)]
s$VAR <- factor(s$VAR, levels=c("YEXO","XPRP","CALO","POPU"))
s$Region <- factor(s$Region, levels=c("WLD","SSA","IND","CHN","SEA","OAS","FSU","MEN","OAM","EUR","ANZ","NAM"))
s <- regionagmiplabelshort(s)
s <- varlabel(s)

s$RCP <- factor(s$RCP, levels=c("RCP85","RCP26"))
s$SSP <- factor(s$SSP, levels=c("SSP2","SSP3"))
interact<-interaction(s$SSP, s$RCP, sep="_")

mins <- ddply(s,c("Region","interact","VAR"),summarize,min=min(Value))
maxs <- ddply(s,c("Region","interact","VAR"),summarize,max=max(Value))
v99 <- ddply(s,c("Region","interact","VAR"),summarize,v99=quantile(Value,0.99))
v01 <- ddply(s,c("Region","interact","VAR"),summarize,v01=quantile(Value,0.01)) 
v05 <- ddply(s,c("Region","interact","VAR"),summarize,v50=mean(Value))
stats <- cbind(mins,maxs,v01,v99,v05)
stats2 <- stats[c(1,2,3,4,8,12,16,20)]

nocc2 <- pdata_long[pdata_long$Year==yearhist & pdata_long$GCM=="nocc", c(-1,-4,-5,-6)]
#nocc2 <- nocc2[(nocc2$VAR=="POPU" | nocc2$VAR=="CALO" | nocc2$VAR=="XPRP"),c(1,2,3)]
nocc2 <- nocc2[(nocc2$Crop=="AGR"& nocc2$VAR=="POPU") | (nocc2$Crop=="AGR"& nocc2$VAR=="CALO") | (nocc2$Crop=="CR5" & nocc2$VAR=="YEXO") | (nocc2$Crop=="AGR"& nocc2$VAR=="XPRP"),c(-3)]
nocc2 <- nocc2[nocc2$Region=="WLD" | nocc2$Region=="CHN" | nocc2$Region=="IND" | nocc2$Region=="FSU" | nocc2$Region=="MEN" | 
                 nocc2$Region=="SSA" | nocc2$Region=="SEA" | nocc2$Region=="OAS" | nocc2$Region=="OAM" |
                 nocc2$Region=="EUR" | nocc2$Region=="ANZ" | nocc2$Region=="NAM" 
               ,c(1,2,3,4)]
nocc2 <- regionagmiplabelshort(nocc2)
nocc2 <- varlabel(nocc2)
nocc2ssp2 <- nocc2[nocc2$SSP=="SSP2",c(-2)]
nocc2ssp3 <- nocc2[nocc2$SSP=="SSP3",c(-2)]

gboxregssp <- ggplot() + 
  geom_boxplot(data=s, aes(x=interact, y=Value, fill=interact)) +
  scale_fill_manual(name="SSP_RCP",values=col3.2) +
  ylab("") +
  xlab("") +	
  MyThemeLine +
  theme(
    axis.text.x = element_text(hjust=1,size = fsize, angle = 90),
    strip.text = element_text(size=fsize, colour = "black", angle = 0,face="plain")
  ) +
  geom_hline(data=nocc2ssp2,aes(yintercept=Value), color="black", size=1) +
  geom_hline(data=nocc2ssp3,aes(yintercept=Value), color="grey70", size=1)

  
  
gboxregssp2<- gboxregssp + facet_grid(VAR~Region,scales = 'free')
plot(gboxregssp2)

filenameboxssp <- paste('../output/png/boxplot2_2050_allreg_4index_N100_noco2SSP.png', sep = '')
ggsave(gboxregssp2,filename=filenameboxssp, dpi=250,width=18, height=12)


#----End --- Boxplot across regions  with co2 and noco2 and RCP2.6 and RCP8.5----

#----- #1-2 Boxplot across region for sensitivity analysis  ----

scensenslabel <- function(y){
  levels(y$SSP)[levels(y$SSP)=="SSP2_TRD2"] <- "SSP2"
  z <- y
  return(z)
}

col6<- c('orange','red','pink','cyan','blue','violet')

yearhist <- c("2050")

s <- pdata[(pdata$SSP=="SSP2_TRD1" | pdata$SSP=="SSP2_TRD2" | pdata$SSP=="SSP2_TRD3") & pdata$Year==yearhist & pdata$CO2fert=="noco2", c(-1,-6)]
s <- s[(s$Crop=="AGR"& s$VAR=="POPU") | (s$Crop=="AGR"& s$VAR=="CALO") | (s$Crop=="CR5" & s$VAR=="YEXO") | (s$Crop=="AGR"& s$VAR=="XPRP"), c(1,2,3,4,5,6,7,8,9)]
s <- s[s$Region=="WLD" | s$Region=="CHN" | s$Region=="IND" | s$Region=="FSU" | s$Region=="MEN" | 
         s$Region=="SSA" | s$Region=="SEA" | s$Region=="OAS" | s$Region=="OAM"
       #         s$Region=="EUR" | s$Region=="ANZ" | s$Region=="NAM" 
       ,c(1,2,3,4,5,6,7,8,9)]
s$VAR <- factor(s$VAR, levels=c("YEXO","XPRP","CALO","POPU"))
s$Region <- factor(s$Region, levels=c("WLD","SSA","IND","CHN","SEA","OAS","FSU","MEN","OAM","EUR","ANZ","NAM"))
s <- regionagmiplabelshort(s)
s <- varlabel(s)
s <- scensenslabel(s)
s$SSP <- factor(s$SSP, levels=c("SSP2_TRD1","SSP2","SSP2_TRD3"))

s$RCP <- factor(s$RCP, levels=c("RCP85","RCP26"))
#s$CO2fert <- factor(s$CO2fert, levels=c("noco2","co2"))
interact<-interaction(s$SSP, s$RCP, sep="_")

mins <- ddply(s,c("Region","interact","VAR"),summarize,min=min(Value))
maxs <- ddply(s,c("Region","interact","VAR"),summarize,max=max(Value))
v99 <- ddply(s,c("Region","interact","VAR"),summarize,v99=quantile(Value,0.99))
v01 <- ddply(s,c("Region","interact","VAR"),summarize,v01=quantile(Value,0.01)) 
v05 <- ddply(s,c("Region","interact","VAR"),summarize,v50=mean(Value))
stats <- cbind(mins,maxs,v01,v99,v05)
stats2 <- stats[c(1,2,3,4,8,12,16,20)]

nocc2 <- pdata_long[pdata_long$Year==yearhist & pdata_long$GCM=="nocc", c(-1,-4,-5,-6)]
#nocc2 <- nocc2[(nocc2$VAR=="POPU" | nocc2$VAR=="CALO" | nocc2$VAR=="XPRP"),c(1,2,3)]
nocc2 <- nocc2[(nocc2$Crop=="AGR"& nocc2$VAR=="POPU") | (nocc2$Crop=="AGR"& nocc2$VAR=="CALO") | (nocc2$Crop=="CR5" & nocc2$VAR=="YEXO") | (nocc2$Crop=="AGR"& nocc2$VAR=="XPRP"),c(-3)]
nocc2 <- nocc2[nocc2$Region=="WLD" | nocc2$Region=="CHN" | nocc2$Region=="IND" | nocc2$Region=="FSU" | nocc2$Region=="MEN" | 
                 nocc2$Region=="SSA" | nocc2$Region=="SEA" | nocc2$Region=="OAS" | nocc2$Region=="OAM"
               #                 nocc2$Region=="EUR" | nocc2$Region=="ANZ" | nocc2$Region=="NAM" 
               ,c(1,2,3,4)]
nocc2 <- regionagmiplabelshort(nocc2)
nocc2 <- varlabel(nocc2)
nocc2ssp2 <- nocc2[nocc2$SSP=="SSP2",c(-2)]
nocc2ssp2trd1 <- nocc2[nocc2$SSP=="SSP2_TRD1",c(-2)]
nocc2ssp2trd3 <- nocc2[nocc2$SSP=="SSP2_TRD3",c(-2)]


gboxregsens <- ggplot() + 
  geom_boxplot(data=s, aes(x=interact, y=Value, fill=interact)) +
  geom_boxplot(data=s, aes(x=interact, y=Value, fill=interact)) +
  scale_fill_manual(name="RCP_CO2fert",values=col6) +
  ylab("") +
  xlab("") +	
  MyThemeLine +
  theme(
    axis.text.x = element_text(hjust=1,size = fsize, angle = 90),
    legend.position = "none",
    strip.text = element_text(size=fsize, colour = "black", angle = 0,face="plain")
  ) +
  geom_hline(data=nocc2ssp2,aes(yintercept=Value), color="black", size=1) +
  geom_hline(data=nocc2ssp2trd1,aes(yintercept=Value), color="grey70", size=0.7) +
  geom_hline(data=nocc2ssp2trd3,aes(yintercept=Value), color="grey70", size=0.7)

gboxregsens2<- gboxregsens + facet_grid(VAR~Region,scales = 'free')
plot(gboxregsens2)


filenameboxsens <- paste('../output/png/boxplot2_2050_region_withhunger_4index_sensrev.png', sep = '')
ggsave(gboxregsens2,filename=filenameboxsens, dpi=250,width=13, height=13)

filenameboxsens <- paste('../output/png/boxplot2_2050_allreg_4index_sens.png', sep = '')
ggsave(gboxregsens2,filename=filenameboxsens, dpi=250,width=18, height=12)


#----End --- Boxplot across regions for sensitivity analysis ----

#-----#2 Ribbon for a single region ----

#---var selection 
varname <- paste("CALO")
paraname <- paste("pcal_fao_tot")
labelname <- paste("Food consumption (kcal/person/day)")
labellog <- paste("Log of food consumption (kcal/person/day)")
cropname <- paste("AGR")

varname <- paste("POPU")
paraname <- paste("risk_fao")
labelname <- paste("Population at risk of hunger (million)")
labellog <- paste("Log of population at risk of hunger (million)")
cropname <- paste("AGR")

#---var selection end


#-- with extreme cases
v <- pdata[pdata$SSP=="SSP2" & pdata$VAR==varname & pdata$Crop==cropname & pdata$CO2fert=="noco2" & pdata$Region=="WLD", c(-2,-8,-9,-10)]

xrib1 <- ddply(v,c("Year","RCP"),transform,max=max(Value))
xrib2 <- ddply(xrib1,c("Year","RCP"),transform,min=min(Value))
xrib3 <- ddply(xrib2,c("Year","RCP"),transform,v1=quantile(Value,0.01))
xrib4 <- ddply(xrib3,c("Year","RCP"),transform,v5=quantile(Value,0.05))
xrib5 <- ddply(xrib4,c("Year","RCP"),transform,v10=quantile(Value,0.1))
xrib6 <- ddply(xrib5,c("Year","RCP"),transform,v35=quantile(Value,0.35))
xrib7 <- ddply(xrib6,c("Year","RCP"),transform,v50=quantile(Value,0.5))
xrib8 <- ddply(xrib7,c("Year","RCP"),transform,v65=quantile(Value,0.65))
xrib9 <- ddply(xrib8,c("Year","RCP"),transform,v90=quantile(Value,0.9))
xrib10 <- ddply(xrib9,c("Year","RCP"),transform,v95=quantile(Value,0.95))
xrib11 <- ddply(xrib10,c("Year","RCP"),transform,v99=quantile(Value,0.99))


xrib11[1] <- lapply(xrib11[1],as.integer)
xrib11[1] <- (xrib11[1]-18)*10+2000

minv <- min(v[7])
maxv <- max(v[7])


#-- noCC case

nocc <- pdata_long[pdata_long$SSP=="SSP2" & pdata_long$GCM=="nocc" & pdata_long$Region=="WLD" & pdata_long$VAR==varname & pdata_long$Crop==cropname, c(-2,-3,-4,-5,-6,-7,-8)]
nocc[1] <- lapply(nocc[1],as.integer)
nocc[1] <- (nocc[1]-18)*10+2000

#-- fao world
fao <- rgdx.param('../data/Risk_fao.gdx',paraname)
#fao <- rgdx.param('../data/Risk_fao.gdx','pcal_fao_tot')
names(fao) <- c("Reg","Year","Value")
fao[3] <- fao[3]*1000 
faoworld <- fao %>% dplyr::filter(Reg == "WLD")
faoworld[2] <- lapply(faoworld[2],as.integer)
faoworld[2] <- (faoworld[2]+1990)


plotrib <-ggplot() +
  geom_ribbon(data=xrib11,aes(x=Year,ymin=min, ymax=max, fill = RCP), alpha=0.1) +
#  geom_ribbon(data=xrib11,aes(x=Year,ymin=v1, ymax=v99, fill = RCP), alpha=0.1) +
#  geom_ribbon(data=xrib11,aes(x=Year,ymin=v5, ymax=v95, fill = RCP), alpha=0.1) +
#  geom_ribbon(data=xrib11,aes(x=Year,ymin=v10, ymax=v90, fill = RCP), alpha=0.2) +
  geom_ribbon(data=xrib11,aes(x=Year,ymin=v35, ymax=v65, fill = RCP), alpha=0.5) +
  geom_line(data=xrib11,aes(x=Year,y=v50,group = RCP, colour=RCP), size = 1) +
  geom_line(data=nocc,aes(x=Year,y=Value), size = 1) +
  geom_line(data=faoworld,aes(x=Year,y=Value), size = 1) +
  ylab(labelname) +
  ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0.02) +
#  labs(aes(colour = RCP)) + 
  scale_colour_manual(name="RCP",values=col) + 
  scale_fill_manual(name="RCP",values=col) +
#  theme_bw(base_size = 30) + labs(y =labelname, x="") +
  MyThemeLine +
  theme(
        legend.position = "none"
        )  

plotrib2 <- plotrib + facet_wrap( ~ RCP,scales="free")

plot(plotrib)

#filenamerib <- paste('../output/png/ribbon_',varname,'_N100.png', sep = '')
#ggsave(plotrib,filename=filenamerib, dpi=250,width=8, height=6)

vbar <- pdata[pdata$SSP=="SSP2" & pdata$Year=="2050" & pdata$VAR==varname & pdata$Crop==cropname & pdata$Region=="WLD", c(-1,-2,-8,-9,-10)] 
vbar$RCP <- factor(vbar$RCP, levels=c("RCP85","RCP26"))
vbar$CO2fert <- factor(vbar$CO2fert, levels=c("noco2","co2"))

vbar2 <- ddply(vbar,c("RCP","CO2fert"),summarise,max=max(Value))
vbar3 <- ddply(vbar,c("RCP","CO2fert"),summarise,min=min(Value))
vbar4 <- ddply(vbar,c("RCP","CO2fert"),summarise,v35=quantile(Value,0.35))
vbar5 <- ddply(vbar,c("RCP","CO2fert"),summarise,v50=quantile(Value,0.5))
vbar6 <- ddply(vbar,c("RCP","CO2fert"),summarise,v65=quantile(Value,0.65))
vbar7 <- ddply(vbar,c("RCP","CO2fert"),summarise,v1=quantile(Value,0.01))
vbar8 <- ddply(vbar,c("RCP","CO2fert"),summarise,v99=quantile(Value,0.99))



vbar9 <- cbind(vbar2,vbar3,vbar4,vbar5,vbar6,vbar7,vbar8)
vbar10 <- vbar9[c(1,2,3,6,9,12,15,18,21)]
ybar <- vbar10
#xrib3co2 <- ddply(xrib2co2,c("RCP","CO2fert"),transform,v35=quantile(Value,0.35))
#xrib4co2 <- ddply(xrib3co2,c("RCP","CO2fert"),transform,v50=quantile(Value,0.5))
#xrib5co2 <- ddply(xrib4co2,c("RCP","CO2fert"),transform,v65=quantile(Value,0.65))

interact<-interaction(ybar$RCP, ybar$CO2fert, sep="_")

gybar <- ggplot()+
  geom_errorbar(data=ybar, aes(x=as.factor(interact), ymin=min, ymax=max,colour=interact), width=.1,size=2,position=position_dodge(.9)) +
  geom_point(data=ybar, aes(x=as.factor(interact), y=v50, colour=interact), shape=16 ,size=3) +
  geom_point(data=ybar, aes(x=as.factor(interact), y=v35, colour=interact), shape=18 ,size=3) +
  geom_point(data=ybar, aes(x=as.factor(interact), y=v65, colour=interact), shape=18 ,size=3) +
  labs(aes(colour = interact)) +
  ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0.02) +
  ylab("") +
  xlab("") +
  MyThemeLine +
  scale_color_manual("Scenarios",values=col3) +
  theme(
#    legend.position = "none",
     axis.text.x=element_blank(),
     axis.text.y=element_blank(),
    axis.line = element_line(colour = "white"),
    axis.text=element_text(size =30,colour = "white"),
    panel.border = element_blank(),
    axis.ticks=element_blank()
  )

plot(gybar)

filenamerib <- paste('../output/png/ribbon_',varname,'_errorbar.png', sep = '')
multfigure2(plotrib,gybar,filenamerib)

multfigure2 <- function(p1,p2,x){
  grid.newpage()
  l <- grid.layout( nrow=1, ncol=2, heights=unit(c(1),c("null")), widths=unit(c(2, 0.75),c("null","null")))
  png(filename=x, width=1800, height=1300,res=200)
  v <- viewport(layout = l)
  pushViewport(v)
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  popViewport()
  dev.off()  
  return()
}

multfigure2pdf <- function(p1,p2,x){
  grid.newpage()
  pdf(file=x)
  l <- grid.layout( nrow=1, ncol=2, heights=unit(c(0.5),c("null")), widths=unit(c(1, 0.75),c("null","null")))
#  l <- grid.layout( nrow=1, ncol=2, heights=unit(c(1),c("null")), widths=unit(c(1, 0.75),c("null","null")))
#  png(filename=x, width=1800, height=1300,res=200)
  v <- viewport(layout = l)
  pushViewport(v)
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  popViewport()
  dev.off()  
  return()
}

filenamerib <- paste('../output/pdf/Figure1ab_',varname,'.pdf', sep = '')
multfigure2pdf(plotrib,gybar,filenamerib)

#-----#2 Ribbon for a single region with and without CO2 fertilization ----

#-- with extreme cases
rcphist <- c("RCP85")
col2<- c('orange','red')

rcphist <- c("RCP26")
col2<- c('cyan','blue')


vco2 <- pdata[pdata$SSP=="SSP2" & pdata$VAR==varname & pdata$Crop==cropname & pdata$RCP==rcphist & pdata$Region=="WLD", c(-2,-5,-8,-9,-10)]

xrib1co2 <- ddply(vco2,c("Year","CO2fert","RCP"),transform,max=max(Value))
xrib2co2 <- ddply(xrib1co2,c("Year","CO2fert","RCP"),transform,min=min(Value))
xrib3co2 <- ddply(xrib2co2,c("Year","CO2fert","RCP"),transform,v35=quantile(Value,0.35))
xrib4co2 <- ddply(xrib3co2,c("Year","CO2fert","RCP"),transform,v50=quantile(Value,0.5))
xrib5co2 <- ddply(xrib4co2,c("Year","CO2fert","RCP"),transform,v65=quantile(Value,0.65))

xrib5co2[1] <- lapply(xrib5co2[1],as.integer)
xrib5co2[1] <- (xrib5co2[1]-18)*10+2000


plotribco2 <-ggplot() +
  geom_ribbon(data=xrib5co2,aes(x=Year,ymin=min, ymax=max, fill = CO2fert), alpha=0.1) +
  geom_ribbon(data=xrib5co2,aes(x=Year,ymin=v35, ymax=v65, fill = CO2fert), alpha=0.5) +
  geom_line(data=xrib5co2,aes(x=Year,y=v50,group = CO2fert, colour=CO2fert), size = 1) +
  geom_line(data=nocc,aes(x=Year,y=Value), size = 1) +
  geom_line(data=faoworld,aes(x=Year,y=Value), size = 1) +
  ylab(labelname) +
  #  labs(aes(colour = RCP)) + 
  scale_colour_manual(name="CO2fert",values=col2) + 
  scale_fill_manual(name="CO2fert",values=col2) +
  #  theme_bw(base_size = 30) + labs(y =labelname, x="") +
  MyThemeLine

plot(plotribco2)

filenameribco2 <- paste('../output/png/ribbon_',varname,rcphist,'_co2fert_N100.png', sep = '')
ggsave(plotribco2,filename=filenameribco2, dpi=250,width=8, height=6)

#-----#2 Ribbon for a trade elas sensitivity analysis ----


vtrd <- pdata[pdata$VAR==varname & pdata$Crop==cropname & pdata$Region=="WLD" & pdata$CO2fert=="noco2" , c(-2,-6,-9,-10)]
vtrd <- vtrd[vtrd$SSP=="SSP2_TRD1" | vtrd$SSP=="SSP2_TRD2"| vtrd$SSP=="SSP2_TRD3", c(1,2,3,4,5,6,7) ]

xrib1trd <- ddply(vtrd,c("Year","SSP","RCP"),transform,max=max(Value))
xrib2trd <- ddply(xrib1trd,c("Year","SSP","RCP"),transform,min=min(Value))
xrib3trd <- ddply(xrib2trd,c("Year","SSP","RCP"),transform,v35=quantile(Value,0.35))
xrib4trd <- ddply(xrib3trd,c("Year","SSP","RCP"),transform,v50=quantile(Value,0.5))
xrib5trd <- ddply(xrib4trd,c("Year","SSP","RCP"),transform,v65=quantile(Value,0.65))

xrib5trd[1] <- lapply(xrib5trd[1],as.integer)
xrib5trd[1] <- (xrib5trd[1]-18)*10+2000

plotribtrd <-ggplot() +
  geom_ribbon(data=xrib5trd,aes(x=Year,ymin=min, ymax=max, fill = SSP), alpha=0.1) +
  geom_ribbon(data=xrib5trd,aes(x=Year,ymin=v35, ymax=v65, fill = SSP), alpha=0.5) +
  geom_line(data=xrib5trd,aes(x=Year,y=v50,group = SSP, colour=SSP), size = 1) +
  geom_line(data=nocc,aes(x=Year,y=Value), size = 1) +
  geom_line(data=faoworld,aes(x=Year,y=Value), size = 1) +
  ylab(labelname) +
  #  labs(aes(colour = RCP)) + 
  scale_colour_manual(name="SSP",values=col) + 
  scale_fill_manual(name="SSP",values=col) +
  #  theme_bw(base_size = 30) + labs(y =labelname, x="") +
  MyThemeLine


plotribtrd2 <- plotribtrd + facet_wrap( ~ RCP, scales="free")

plot(plotribtrd2)

filenameribtrd <- paste('../output/png/ribbon_',varname,rcphist,'_trd.png', sep = '')
ggsave(plotribtrd,filename=filenameribtrd, dpi=250,width=8, height=6)





#-----#3  Histogram for a single year ----
yearhist <- c("2050")
w <- v[v$Year==yearhist,c(-1)]
nocc1 <- nocc[nocc$Year==yearhist,c(-1)]

minv <- min(min(w[6]),nocc1)
maxv <- max(max(w[6]),nocc1)
muv <- ddply(w,"RCP", summarise, grp.mean=mean(Value))

plothist <- ggplot() + 
  #  geom_histogram(data=w, aes(x=Value,fill=RCP,y=..density..,color=RCP),alpha = 0.3, position = "identity") +
  geom_density(data=w, aes(x=Value,fill=RCP,y=..density..,color=RCP),alpha = 0.3, size=1, position = "identity") +  
  xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
#  xlim(0, 200) +
  scale_fill_manual(name="RCP",values=col) +
  scale_color_manual(name="RCP",values=col) +
  ylab("Density") + 
  xlab(labelname) +
  ggtitle(yearhist) +
  MyThemeLine +
  geom_vline(data=muv,aes(xintercept=grp.mean, color=RCP), linetype="dashed", size=1) +
  geom_vline(aes(xintercept=nocc1), color="black", size=1)

plot(plothist)

filenamehist <- paste('../output/png/histogram_',varname,'_',yearhist,'_noco2_N100.png', sep = '')
ggsave(plothist,filename=filenamehist, dpi=250,width=8, height=6)


plotcdf <- ggplot() + 
  geom_step(data=w, aes(x=Value,y=..y..,color=RCP),stat="ecdf") +
  #  scale_x_log10() +
  scale_fill_manual(name="RCP",values=col) +
  scale_color_manual(name="RCP",values=col) +
  ggtitle(yearhist) +
  ylab("Density") + 
  xlab(labelname) +
  MyThemeLine +
  geom_vline(data=muv,aes(xintercept=grp.mean, color=RCP), linetype="dashed", size=1) +
  geom_vline(aes(xintercept=nocc1), color="black", size=1) +
  geom_hline(yintercept=0.1, color="black", linetype = "dashed") +
  geom_hline(yintercept=0.01, color="black", linetype = "dashed")

plot(plotcdf)
filenamecdf <- paste('../output/png/CDF_',varname,'_',yearhist,'_noco2_N100.png', sep = '')
ggsave(plotcdf,filename=filenamecdf, dpi=250,width=8.5, height=6)


#-----#3  Histogram for a single year with RCO26 and RCP85 and co2 and nonco2 ----
yearhist <- c("2050")
w <- pdata[pdata$SSP=="SSP2" & pdata$Year==yearhist & pdata$VAR==varname & pdata$Crop==cropname & pdata$Region=="WLD", c(-1,-2,-8,-9,-10)]
nocc1 <- nocc[nocc$Year==yearhist,c(-1)]

w$RCP <- factor(w$RCP, levels=c("RCP85","RCP26"))
w$CO2fert <- factor(w$CO2fert, levels=c("noco2","co2"))
interact<-interaction(w$RCP, w$CO2fert, sep="_")

minv <- min(min(w[6]),nocc1)
maxv <- max(max(w[6]),nocc1)
muv <- ddply(w,"interact", summarise, grp.mean=mean(Value))

plothistrcpco2 <- ggplot() + 
  #  geom_histogram(data=w, aes(x=Value,fill=RCP,y=..density..,color=RCP),alpha = 0.3, position = "identity") +
  geom_density(data=w, aes(x=Value,fill=interact,y=..density..,color=interact),alpha = 0.3, size=1, position = "identity") +  
  xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
  xlim(300, 800) +
  scale_fill_manual(name="Scenarios",values=col3) +
  scale_color_manual(name="Scenarios",values=col3) +
  ylab("Density") + 
  xlab(labelname) +
  ggtitle(yearhist) +
  MyThemeLine +
  geom_vline(data=muv,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1) +
  geom_vline(aes(xintercept=nocc1), color="black", size=1)

plot(plothistrcpco2)

filenamehistrcpco2 <- paste('../output/png/histogram_',varname,'_',yearhist,'_N100_rcpco2s_legend.pdf', sep = '')
ggsave(plothistrcpco2,filename=filenamehistrcpco2, dpi=250,width=8, height=6)

filenamehistrcpco2 <- paste('../output/pdf/Figure1cd_',varname,'_',yearhist,'.pdf', sep = '')
ggsave(plothistrcpco2,filename=filenamehistrcpco2, dpi=250,width=8, height=6)


#-----#3  Histogram for a single year across regions with RCO26 and RCP85 and co2 and nonco2----

#Making multi figures
multifigure11 <- function(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,x){
  grid.newpage()
  png(filename=x, width=2000, height=2500,res=200)
  l <- grid.layout(nrow=4, ncol=3, heights=unit(c(1,1,1,1),c("null","null","null","null")), widths=unit(c(1),c("null")),respect = FALSE)
  v <- viewport(layout = l)
  pushViewport(v)
  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=1,layout.pos.col=3))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p5, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  print(p6, vp = viewport(layout.pos.row=2,layout.pos.col=3))
  print(p7, vp = viewport(layout.pos.row=3,layout.pos.col=1))
  print(p8, vp = viewport(layout.pos.row=3,layout.pos.col=2))
  print(p9, vp = viewport(layout.pos.row=3,layout.pos.col=3))
  print(p10, vp = viewport(layout.pos.row=4,layout.pos.col=1))
  print(p11, vp = viewport(layout.pos.row=4,layout.pos.col=2))
  
  popViewport()
  dev.off()  
  return()
}

multifigure8 <- function(p1,p2,p3,p4,p5,p6,p7,p8,x){
  grid.newpage()
  png(filename=x, width=2000, height=2000,res=200)
  l <- grid.layout(nrow=3, ncol=3, heights=unit(c(1,1,1),c("null","null","null")), widths=unit(c(1),c("null")),respect = FALSE)
  v <- viewport(layout = l)
  pushViewport(v)
  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=1,layout.pos.col=3))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p5, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  print(p6, vp = viewport(layout.pos.row=2,layout.pos.col=3))
  print(p7, vp = viewport(layout.pos.row=3,layout.pos.col=1))
  print(p8, vp = viewport(layout.pos.row=3,layout.pos.col=2))
  popViewport()
  dev.off()  
  return()
}


varname <- paste("YEXO")
labelname <- paste("yield (2005=1)")
cropname <- paste("CR5")

varname <- paste("XPRP")
labelname <- paste("Price index (baseyear=1)")
cropname <- paste("AGR")

varname <- paste("CALO")
labelname <- paste("Food consumption\n(kcal/person/day)")
cropname <- paste("AGR")

varname <- paste("POPU")
labelname <- paste("Population at risk of\nhunger (million)")
cropname <- paste("AGR")

yearhist <- c("2050")
region_list <- c("WLD","CHN","IND","FSU","MEN","SSA","SEA","OAS","OAM","EUR","ANZ","NAM")
region_list <- c("WLD","CHN","IND","FSU","MEN","SSA","SEA","OAS","OAM")
#region_list <- c("SSA")
#region_list <- c("CHN")


w0 <- pdata[pdata$SSP=="SSP2" & pdata$Year==yearhist & pdata$VAR==varname & pdata$Crop==cropname, c(-1,-8,-9,-10)]
nocc0 <- pdata_long[pdata_long$Year==yearhist & pdata_long$SSP=="SSP2" & pdata_long$GCM=="nocc" & pdata_long$VAR==varname & pdata_long$Crop==cropname, c(-1,-3,-4,-5,-6,-7,-8)]
minv <- min(min(w0[7]),min(nocc0[2]))
maxv <- max(max(w0[7]),max(nocc0[2]))

for(r in 1:length(region_list)){
  nocc1 <- 0
  regionname <- region_list[r]
  regionfullname <- fregionagmipname(regionname)

  w <- pdata[pdata$SSP=="SSP2" & pdata$Year==yearhist & pdata$VAR==varname & pdata$Crop==cropname & pdata$Region==regionname, c(-1,-2,-8,-9,-10)]
  nocc1 <- pdata_long[pdata_long$Year==yearhist & pdata_long$SSP=="SSP2" & pdata_long$GCM=="nocc" & pdata_long$VAR==varname & pdata_long$Crop==cropname & pdata_long$Region==regionname, c(-1,-2,-3,-4,-5,-6,-7,-8)]
  
  w$RCP <- factor(w$RCP, levels=c("RCP85","RCP26"))
  w$CO2fert <- factor(w$CO2fert, levels=c("noco2","co2"))
  interact<-interaction(w$RCP, w$CO2fert, sep="_")
  
  muv <- ddply(w,"interact", summarise, grp.mean=mean(Value))
  minv <- min(min(w[6]),nocc1)
  maxv <- max(max(w[6]),nocc1)
  
  plothist <- ggplot() + 
  #  geom_histogram(data=w, aes(x=Value,fill=RCP,y=..density..,color=RCP),alpha = 0.3, position = "identity") +
    geom_density(data=w, aes(x=Value,fill=interact,y=..density..,color=interact),alpha = 0.3, size=1, position = "identity") +  
    xlim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) +
#    xlim(2500, 4000) +
#    xlim(0, 250) +
    scale_fill_manual(name="Scenarios",values=col3) +
    scale_color_manual(name="Scenarios",values=col3) +
    ylab("Density") + 
    xlab(labelname) +
    ggtitle(regionfullname) +
    MyThemeLine +
    theme(legend.position = "none") +
    geom_vline(data=muv,aes(xintercept=grp.mean, color=interact), linetype="dashed", size=1) +
    annotate("segment",xend=nocc1,x=nocc1,y=-Inf,yend=+Inf, color="black", size=1)
#    geom_vline(aes(xintercept=nocc1), color="black", size=1)

  plot(plothist)

  if (regionname=="CHN"){
    g1 <- plothist
  } else if (regionname=="IND"){
    g2 <- plothist  
  } else if (regionname=="FSU"){
    g3 <- plothist  
  } else if (regionname=="MEN"){
    g4 <- plothist   
  } else if (regionname=="SSA"){
    g5 <- plothist  
  } else if (regionname=="SEA"){
    g6 <- plothist  
  } else if (regionname=="OAS"){
    g7 <- plothist  
  } else if (regionname=="OAM"){
    g8 <- plothist  
  } else if (regionname=="EUR"){
    g9 <- plothist  
  } else if (regionname=="ANZ"){
    g10 <- plothist  
  } else if (regionname=="NAM"){
    g11 <- plothist  
  } else {
  }
  ;
  
}

filenamehist <- paste('../output/png/histogram/',varname,'_',cropname,'_',yearhist,'_reg_rcpco2rev.png', sep = '')
#ggsave(plothist,filename=filenamehist, dpi=250,width=8, height=6)
multifigure11(g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,filenamehist)
multifigure8(g1,g2,g3,g4,g5,g6,g7,g8,filenamehist)

#-----#4  Histogram for a single RCP across years----

rcphist <- c("RCP85")
rcphist <- c("RCP26")
u <- v[v$RCP==rcphist & v$Year!="2010",c(-4)]
#nocc1 <- nocc[nocc$Year==yearhist,c(-1)]

minu <- min(min(u[6]))
maxu <- max(max(u[6]))

muu <- ddply(u,"Year", summarise, grp.mean=mean(Value))

plothistrcp <- ggplot() + 
  #  geom_histogram(data=w, aes(x=Value,fill=RCP,y=..density..,color=RCP),alpha = 0.3, position = "identity") +
  geom_density(data=u, aes(x=Value,fill=Year,y=..density..,color=Year),alpha = 0.3, size=1, position = "identity") +  
  xlim(minu-abs(minu)*0.01, maxu+abs(maxu)*0.01) +
  #  xlim(200, 800) +
  ylim(0, 0.015) +
#  scale_fill_manual(name="Year",values=col) +
#  scale_color_manual(name="Year",values=col) +
  ylab("Density") + 
  xlab(labelname) +
  ggtitle(rcphist) +
  MyThemeLine +
  geom_vline(data=muu,aes(xintercept=grp.mean, color=Year), linetype="dashed", size=1)
#  geom_vline(aes(xintercept=nocc1), color="black", size=1)

plot(plothistrcp)

filenamehistrcp <- paste('../output/png/histogram_',varname,'_',rcphist,'_noco2_N100.pdf', sep = '')
ggsave(plothistrcp,filename=filenamehistrcp, dpi=250,width=8, height=6)
filenamehistrcp <- paste('../output/pdf/Figure2ab_',varname,'_',rcphist,'.pdf', sep = '')
ggsave(plothistrcp,filename=filenamehistrcp, dpi=250,width=8, height=6)


plotcdfrcp <- ggplot() + 
  geom_step(data=u, aes(x=Value,y=..y..,color=Year),stat="ecdf") +
  #  scale_x_log10() +
  scale_fill_manual(name="Year",values=col) +
  scale_color_manual(name="Year",values=col) +
  ggtitle(rcphist) +
  ylab("Density") + 
  xlab(labelname) +
  MyThemeLine +
  geom_vline(data=muu,aes(xintercept=grp.mean, color=Year), linetype="dashed", size=1) +
#  geom_vline(aes(xintercept=nocc1), color="black", size=1) +
  geom_hline(yintercept=0.1, color="black", linetype = "dashed") +
  geom_hline(yintercept=0.01, color="black", linetype = "dashed")

plot(plotcdfrcp)
filenamecdfrcp <- paste('../output/png/CDF_',varname,'_',rcphist,'_noco2_N100.png', sep = '')
ggsave(plotcdfrcp,filename=filenamecdfrcp, dpi=250,width=8.5, height=6)

#----End --- PPOP_UNSH & PCAL: Histogram for a single region----

#-----#4  Histogram for a single RCP across GCMs----
yearhist <- c("2050")
rcphist <- c("RCP26")
rcphist <- c("RCP85")
o <- v[v$RCP==rcphist & v$Year==yearhist,c(-1,-4)]
#nocc1 <- nocc[nocc$Year==yearhist,c(-1)]

mino <- min(min(o[5]))
maxo <- max(max(o[5]))

muo <- ddply(o,"GCM", summarise, grp.mean=mean(Value))

plothistgcm <- ggplot() + 
  #  geom_histogram(data=w, aes(x=Value,fill=RCP,y=..density..,color=RCP),alpha = 0.3, position = "identity") +
  geom_density(data=o, aes(x=Value,fill=GCM,y=..density..,color=GCM),alpha = 0.1, size=1, position = "identity") +  
  xlim(mino-abs(mino)*0.01, maxo+abs(maxo)*0.01) +
#  xlim(200, 800) +
  ylim(0, 0.015) +
  #  scale_fill_manual(name="Year",values=col) +
  #  scale_color_manual(name="Year",values=col) +
  ylab("Density") + 
  xlab(labelname) +
  ggtitle(rcphist, yearhist) +
  MyThemeLine
#  geom_vline(data=muo,aes(xintercept=grp.mean, color=GCM), linetype="dashed", size=1)
#  geom_vline(aes(xintercept=nocc1), color="black", size=1)

plot(plothistgcm)

filenamehistgcm <- paste('../output/png/histogram_',varname,'_',rcphist,'_',yearhist,'_noco2_gcm_N100.png', sep = '')
ggsave(plothistgcm,filename=filenamehistgcm, dpi=250,width=8, height=6)


plotcdfgcm <- ggplot() +
  geom_step(data=o, aes(x=Value,y=..y..,color=GCM),stat="ecdf") +
  #  scale_x_log10() +
#  scale_fill_manual(name="GCM",values=col) +
#  scale_color_manual(name="GCM",values=col) +
  ggtitle(rcphist,yearhist) +
  ylab("Density") + 
  xlab(labelname) +
  MyThemeLine
#  geom_vline(data=muo,aes(xintercept=grp.mean, color=GCM), linetype="dashed", size=1) +
#  geom_vline(aes(xintercept=nocc1), color="black", size=1) +
#  geom_hline(yintercept=0.1, color="black", linetype = "dashed") +
#  geom_hline(yintercept=0.01, color="black", linetype = "dashed")

plot(plotcdfgcm)
filenamecdfgcm <- paste('../output/png/CDF_',varname,'_',rcphist,'_',yearhist,'_noco2_gcm_N100.png', sep = '')
ggsave(plotcdfgcm,filename=filenamecdfgcm, dpi=250,width=8.5, height=6)

#----End --- PPOP_UNSH & PCAL: Histogram for a single region----


#----#3 Food stock and return period
var_list <- c("CALO","POPU")
ssp_list <- c("SSP2","SSP3")
rcp_list <- c("RCP85","RCP26")
CO2fert_list <- c("noco2")

region_list <- c("WLD","CAN","USA","BRA","CHN","IND","OSA","FSU","EUR","MEN","SSA","SEA","OAS","ANZ","NAM","OAM","AME","SAS")
#region_list <- c("WLD","CHN","IND","FSU","EUR","MEN","SSA","SEA","OAS","ANZ","NAM","OAM")

year_list <- c("2020","2030","2040","2050")
#year_list <- c("2050")


#names(z) <- c("VAR","Year","Region","SSP","RCP","CO2fert","Value")	#Rename the column

yall <- 0
yall <- data.frame(1,1,1,1,1,1,1,1)	#Rename the column
names(yall) <- c("VAR","Year","Region","RCP","SSP","CO2fert","V","return")	#Rename the column

returnp<-c(0,0.01,0.05,0.1,0.35,0.5,0.65,0.9,0.95,0.99,1)
returnn<-c("MIN","RP001","RP005","RP01","RP035","RP05","RP065","RP09","RP095","RP099","MAX")
for(v in 1:length(var_list)){
  for(yr in 1:length(year_list)){
    for(ssp in 1 :length(ssp_list)){
      for (rc in 1 : length(rcp_list)){
        for(r in 1:length(region_list)){
          for(co2 in 1:length(CO2fert_list)){
          for(i in 1 : 11){
            
            varname <- var_list[v]
            yearname <- year_list[yr]
            sspname <- ssp_list[ssp]
            rcpname <- rcp_list[rc]
            regionname <- region_list[r]
            co2fertname <- CO2fert_list[co2]
            cropname <- c("AGR")

            fl <- pdata[pdata$VAR==varname & pdata$Crop==cropname & pdata$CO2fert==co2fertname & pdata$SSP==sspname & pdata$RCP==rcpname & pdata$Year==yearname & pdata$Region==regionname, c(11)]
          
            if (length(fl)>0){
              a <- returnp[i]
              b <- returnn[i]
              v1 <- quantile(fl,a)
              y5 <- data.frame(varname,year_list[yr],region_list[r],rcp_list[rc],ssp_list[ssp],CO2fert_list[co2],b,v1[1])	#Rename the column
              names(y5) <- c("VAR","Year","Region","RCP","SSP","CO2fert","V","return")	#Rename the column
                yall <- rbind(yall,y5) 
            y5<-0
            }
            }
      }
    }
  }
 }
}
}

symDim <- 8
attr(yall, "symName") <- "Foodshock"
lst <- wgdx.reshape(yall, symDim)
#wgdx.reshape(yall, symDim, gdxName = "../data/FoodshockR.gdx")
#wgdx.reshape(yall, symDim, gdxName = "../data/FoodshockR_2050noco2.gdx")
#wgdx.reshape(yall, symDim, gdxName = "../data/FoodshockR_2050co2.gdx")
system(paste("gams","../prog/foodstock.gms",sep=" "))

ReturnPeriodFS <- rgdx.param('../output/gdx/Pfoodstock.gdx','Pfoodstock_agg')
ReturnPeriodpr <- rgdx.param('../output/gdx/Pfoodstock.gdx','PfoodstockRP_agg')
ReturnPeriodVal <- rgdx.param('../output/gdx/Pfoodstock.gdx','PfoodstockVal_agg')
ReturnPeriodValpr <- rgdx.param('../output/gdx/Pfoodstock.gdx','PfoodstockRPVal_agg')

ReturnPeriodpr <- ReturnPeriodpr[c(7)]
ReturnPeriodValpr <- ReturnPeriodValpr[c(7)]
#ReturnPeriodFS[7] <- ReturnPeriodFS[7]*100
ReturnPeriod <- cbind(ReturnPeriodFS,ReturnPeriodpr)
ReturnPeriodVal <- cbind(ReturnPeriodVal,ReturnPeriodValpr)

colnamesRP0 <- c("Year","Region","RCP","SSP","CO2fert","RetPerI","Stock","RetPer")
names(ReturnPeriod) <- colnamesRP0
names(ReturnPeriodVal) <- colnamesRP0

actualstock <- rgdx.param('../output/gdx/Pfoodstock.gdx','ActualStock')
colnamesAS <- c("Region","Stock")
names(actualstock) <- colnamesAS

sspname <- paste("SSP2")
#sspname <- paste("SSP3")
yearname <- paste("2050")

colnamesRP <- c("Year","Region","RCP","SSP","CO2fert","RetPerI","Stock","RetPer")

x0 <- ReturnPeriod[ReturnPeriod$Year==yearname & ReturnPeriod$RCP!="nocc" & ReturnPeriod$SSP==sspname,c(-1,-4,-6)]

names(x0) <- c("Region","RCP","CO2fert","Stock","RetPer")
x <- regionlabelagg(x0)
x$RCP <- factor(x$RCP, levels=c("RCP85","RCP26"))
x$CO2fert <- factor(x$CO2fert, levels=c("noco2","co2"))

xnwld=x[x$Region!="World",c(1,2,3,4,5)]
interact<-interaction(xnwld$RCP, xnwld$CO2fert, sep="_")

plot1 <- ggplot() + 
  geom_line(data=xnwld, aes(x=RetPer,y=Stock,group=interact,color=interact)) + MyThemeLine +
  scale_color_manual(name="interact",values=col3) +
  xlab("Return period (Year)") +
  ylab("Food requirement (million tonne cereal equivalent)")
plot2 <- plot1+ facet_wrap( ~ Region,ncol=3,scales="free")
plot(plot2)
filename <- paste("../output/png/return_stock_region_N100rev.png")
ggsave(plot2,filename=filename, dpi=250,width=10, height=7)

xwld=x[x$Region=="World",c(1,2,3,4,5)]
interact<-interaction(xwld$RCP, xwld$CO2fert, sep="_")

plotwld <- ggplot() + 
  geom_line(data=xwld, aes(x=RetPer,y=Stock,group=interact,colour=interact)) + MyThemeLine +
  scale_color_manual(name="Scenarios",values=col3) +
  xlab("Return period (Year)") +
  ylab("Food requirement \n(million tonne cereal equivalent)") +
  theme(legend.position = c(0.7,0.22))

plot(plotwld)
filenamewld <- paste("../output/png/return_stock_world_N100rev.pdf")
ggsave(plotwld,filename=filenamewld, dpi=250,width=5, height=5)
filenamewld <- paste("../output/pdf/Figure4a.pdf")
ggsave(plotwld,filename=filenamewld, dpi=250,width=5, height=5)


#----#3 Food stock and return period SSP2 and SSP3

colnamesRP <- c("Year","Region","RCP","SSP","CO2fert","RetPerI","Stock","RetPer")

xssp0 <- ReturnPeriod[ReturnPeriod$Year==yearname & ReturnPeriod$CO2fert=="noco2" & ReturnPeriod$RCP!="nocc",c(-1,-5,-6)]

names(xssp0) <- c("Region","RCP","SSP","Stock","RetPer")
xssp <- regionlabelagg(xssp0)

xssp$SSP <- factor(xssp$SSP, levels=c("SSP2","SSP3"))
xssp$RCP <- factor(xssp$RCP, levels=c("RCP85","RCP26"))


xsspnwld=xssp[xssp$Region!="World",c(1,2,3,4,5)]
interact<-interaction(xsspnwld$SSP, xsspnwld$RCP, sep="_")

plotssp1 <- ggplot() + 
  geom_line(data=xsspnwld, aes(x=RetPer,y=Stock,group=interact,color=interact)) + MyThemeLine +
  scale_color_manual(name="interact",values=col3) +
  xlab("Return period (Year)") +
  ylab("Food requirement (million tonne cereal equivalent)")
plotssp2 <- plotssp1+ facet_wrap( ~ Region,ncol=3,scales="free")
plot(plotssp2)

filename <- paste("../output/png/return_stock_region_ssp_noco2rev.png")
ggsave(plotssp2,filename=filename, dpi=250,width=10, height=7)

xsspwld=xssp[xssp$Region=="World",c(1,2,3,4,5)]
interact<-interaction(xsspwld$SSP, xsspwld$RCP, sep="_")

plotwldssp <- ggplot() + 
  geom_line(data=xsspwld, aes(x=RetPer,y=Stock,group=interact,colour=interact)) + MyThemeLine +
  scale_color_manual(name="Scenarios",values=col3) +
  xlab("Return period (Year)") +
  ylab("Food requirement\n(million tonne cereal equivalent)") +
  theme(legend.position = c(0.7,0.22))

plot(plotwldssp)

filename <- paste("../output/png/return_stock_world_ssp_noco2rev.png")
ggsave(plotwldssp,filename=filename, dpi=250,width=5, height=5)

#----stacked Stock vs return period

x1=x[x$Region!="World" & x$RCP=="RCP85" & x$CO2fert=="noco2",c(-2,-3)]

plot3 <- ggplot() + 
  geom_area(data=x1, aes(x=RetPer, y=Stock, fill = Region), stat = "identity") +   MyThemeLine + 
  scale_fill_manual(values =linepalette) +
  xlab("Return period (Year)")+
  ylab("Food requirement\n(million tonne cereal equivalent)") + 
  theme(legend.position = "none")
plot(plot3)

filename <- paste("../output/png/stacked_return_stock_reg_N100rev.pdf")
ggsave(plot3,filename=filename, dpi=250,width=10, height=6)
filename <- paste("../output/pdf/Figure4b.pdf")
ggsave(plot3,filename=filename, dpi=250,width=10, height=6)

#----stacked Stock vs return period for monetary value of stocks

xval0 <- ReturnPeriodVal[ReturnPeriodVal$Year==yearname & ReturnPeriodVal$RCP!="nocc" & ReturnPeriodVal$SSP==sspname,c(-1,-4,-6)]

names(xval0) <- c("Region","RCP","CO2fert","Stock","RetPer")
xval <- regionlabelagg(xval0)
xval$RCP <- factor(xval$RCP, levels=c("RCP85","RCP26"))
xval$CO2fert <- factor(xval$CO2fert, levels=c("noco2","co2"))

xval1=xval[xval$Region!="World" & xval$RCP=="RCP85" & xval$CO2fert=="noco2",c(-2,-3)]

plotval3 <- ggplot() + 
  geom_area(data=xval1, aes(x=RetPer, y=Stock, fill = Region), stat = "identity") +   MyThemeLine + 
  scale_fill_manual(values =linepalette) +
  xlab("Return period (Year)")+
  ylab("Food requirement\n(million USD)") + 
  theme(legend.position = "none")
plot(plotval3)

filename <- paste("../output/png/stacked_return_stock_reg_N100_valrev.pdf")
ggsave(plotval3,filename=filename, dpi=250,width=10, height=6)
filename <- paste("../output/pdf/Figure4d.pdf")
ggsave(plotval3,filename=filename, dpi=250,width=10, height=6)

#--- Chart comparing actal & required stocks
reqstock <- x0[x0$Region!="WLD" & x0$RCP=="RCP85" & x0$RetPer==100 & x0$CO2fert=="noco2",c(1,4)]
reqstock[3] <- c("Requirement")
actualstock1 <- actualstock[actualstock$Region!="WLD",c(1,2)]
actualstock1[3] <- c("Actual \n reverse")

stocktable <- rbind(reqstock,actualstock1)
stocktable1 <- regionlabelagg(stocktable) 

gbar_stock <- ggplot() + 
  geom_bar(data=stocktable1, aes(x=as.factor(V3),y=Stock,fill=Region), position="stack", stat="identity", colour='black') +
  ylab("Food requirement\n(million tonne cereal equivalent)") + xlab("") +	MyThemeLine +
  scale_fill_manual(values =linepalette) + theme(axis.text.x=element_text(hjust=0.5))
  
plot(gbar_stock)  
filenamestack <- paste("../output/png/stached_bar_stockcomparison_N100rev.pdf")
ggsave(gbar_stock,filename=filenamestack, dpi=250,width=7, height=6)
filenamestack <- paste("../output/pdf/Figure4c.pdf")
ggsave(gbar_stock,filename=filenamestack, dpi=250,width=7, height=6)


#--- Combine food stock figures

#filenamefoodstock <- paste("../output/png/food_stock_world_N100.png")
filenamefoodstock <- paste("../output/png/food_stock_world_Valaddedrev2.png")

multifigure4 <- function(p1,p2,p3,p4,x){
  grid.newpage()
  png(filename=x, width=2000, height=1600,res=200)
  l <- grid.layout(nrow=2, ncol=2, heights=unit(c(1,1),c("null","null")), widths=unit(c(1,1),c("null","null")),respect = FALSE)
  v <- viewport(layout = l)
  pushViewport(v)
#  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  
  popViewport()
  dev.off()  
  return()
}

multifigure4(plotwld,gbar_stock,plot3,plotval3,filenamefoodstock)

multifigure4pdf <- function(p1,p2,p3,p4,x){
  grid.newpage()
  pdf(file=x)
#  png(filename=x, width=2000, height=1600,res=200)
  l <- grid.layout(nrow=2, ncol=2, heights=unit(c(1,1),c("null","null")), widths=unit(c(2,2),c("null","null")),respect = FALSE)
  v <- viewport(layout = l)
  pushViewport(v)
  #  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=2,layout.pos.col=1))
  print(p4, vp = viewport(layout.pos.row=2,layout.pos.col=2))
  
  popViewport()
  dev.off()  
  return()
}

multifigure4pdf(plotwld,gbar_stock,plot3,plotval3,filenamefoodstock)


multifigure42 <- function(p1,p2,p3,p4,x){
  grid.newpage()
  png(filename=x, width=2600, height=800,res=200)
  l <- grid.layout(nrow=1, ncol=4, heights=unit(c(1),c("null")), widths=unit(c(1,0.8,0.8,1),c("null","null","null","null")),respect = FALSE)
  v <- viewport(layout = l)
  pushViewport(v)
  grid.text(x2, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2))
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  print(p3, vp = viewport(layout.pos.row=1,layout.pos.col=3))
  print(p4, vp = viewport(layout.pos.row=1,layout.pos.col=4))
  
  popViewport()
  dev.off()  
  return()
}

multifigure42(plotwld,plot3,plotval3,gbar_stock,filenamefoodstock)

#---Figure1-b/c
regionname <- paste("WLD")

x3 <- ReturnPeriod[ReturnPeriod$Year=="2050" & ReturnPeriod$Region==regionname &  ReturnPeriod$CO2fert=="noco2",c(-1,-2,-4,-5,-6)]
names(x3) <- c("RCP","Stock","RetPer")
colrp <- c('red','blue')
plot1 <- ggplot() + 
  geom_line(data=x3, aes(x=RetPer,y=Stock,group=RCP,color=RCP), size=1.5) + MyThemeLine +
  scale_color_manual(values =colrp) +
  xlab("Return period (Year)") +
  ylab(ylabelname)
plot(plot1)
filename <- paste("../output/png/stackedd_return_stock_WLD_N100.png")
ggsave(plot1,filename=filename, dpi=250,width=8.5, height=6)

#Figure1-c/b SDM

varname <- paste("POPU")

#!!!!! SELECT !!!!!
#SDM0 <- rgdx.param('../output/gdx/Pfoodstock.gdx','PSDM_nocc')
SDM0 <- rgdx.param('../output/gdx/Pfoodstock.gdx','PSDM')
#!!!!! SELECT !!!!!
colnamesRP0 <- c("VAR","Year","Region","RCP","SSP","CO2fert","RetPerI","Stock")
names(SDM0) <- colnamesRP0
SDM  <- SDM0[SDM0$VAR==varname, c(-1)]
SDM$RCP <- factor(SDM$RCP, levels=c("RCP85","RCP26"))
xsdm <- SDM[(SDM$RetPerI=="RP05" | SDM$RetPerI=="RP065" | SDM$RetPerI=="RP09" | SDM$RetPerI=="RP095" | SDM$RetPerI=="RP099"| SDM$RetPerI=="MAX" ) & SDM$Region=="WLD" &  SDM$CO2fert=="noco2"&  SDM$SSP=="SSP2" ,c(-2,-4,-5,-8)]
#names(xsdm) <- c("Year","RCP","Stock")

interact<-interaction(xsdm$RCP, xsdm$RetPerI, sep=" : ")
colsdm<- c('lightblue','dodgerblue','blue2','blue3','blue4','black')

RPlabel <- function(y){
  levels(y$RetPerI)[levels(y$RetPerI)=="RP05"]  <- "2"
  levels(y$RetPerI)[levels(y$RetPerI)=="RP065"] <- "3"
  levels(y$RetPerI)[levels(y$RetPerI)=="RP09"]  <- "10"
  levels(y$RetPerI)[levels(y$RetPerI)=="RP095"] <- "20"
  levels(y$RetPerI)[levels(y$RetPerI)=="RP099"] <- "100"
  levels(y$RetPerI)[levels(y$RetPerI)=="MAX"] <- ">100"
  z <- y
  return(z)
}

xsdm <- RPlabel(xsdm)
xsdm$RetPerI <- factor(xsdm$RetPerI, levels=c("2","3","10","20","100",">100"))

plot2 <- ggplot() + 
  geom_line(data=xsdm, aes(x=Year,y=Stock, group=interact, color=RetPerI), size=1.5) +  MyThemeLine +
  scale_color_manual("Return period",values =colsdm) +
  xlab("") + 
#!!!!! SELECT !!!!!
#  ylab("Population at risk of hunger\n [% change from the no climate change level]")
  ylab("Population at risk of hunger [% change\nfrom the median climate state]")
#!!!!! SELECT !!!!!

plot3 <- plot2 + facet_wrap( ~ RCP) 
plot(plot3)

#!!!!! SELECT !!!!!
#filenamesdm <- paste('../output/png/return_SDM_vsNoCC.png')
filenamesdm <- paste('../output/png/return_SDM_vsMean_N100.pdf')
#!!!!! SELECT !!!!!
ggsave(plot3,filename=filenamesdm, dpi=250,width=7.5, height=3.5)

filenamesdm <- paste('../output/pdf/Figure2cd.pdf')
ggsave(plot3,filename=filenamesdm, dpi=250,width=7.5, height=3.5)

#-- by regions

xsdmreg <- SDM[(SDM$RetPerI=="RP05" | SDM$RetPerI=="RP065" | SDM$RetPerI=="RP09" | SDM$RetPerI=="RP095" | SDM$RetPerI=="RP099"| SDM$RetPerI=="MAX" ) &  SDM$CO2fert=="noco2" ,c(-4,-5,-8)]
xsdmreg <- xsdmreg[xsdmreg$Region=="CHN" | xsdmreg$Region=="IND" | xsdmreg$Region=="FSU" | xsdmreg$Region=="MEN" | 
                     xsdmreg$Region=="SSA" | xsdmreg$Region=="SEA" | xsdmreg$Region=="OAS" | xsdmreg$Region=="OAM" |
                     xsdmreg$Region=="EUR" | xsdmreg$Region=="ANZ" | xsdmreg$Region=="NAM" ,c(1,2,3,4,5)]

interact<-interaction(xsdmreg$RCP, xsdmreg$RetPerI, sep=" : ")
xsdmreg <- RPlabel(xsdmreg)
xsdmreg <- regionagmiplabelshort(xsdmreg)
xsdmreg$RetPerI <- factor(xsdmreg$RetPerI, levels=c("2","3","10","20","100","MAX"))

plot4 <- ggplot() + 
  geom_line(data=xsdmreg, aes(x=Year,y=Stock, group=interact, color=RetPerI), size=1.5) + MyThemeLine +
  scale_color_manual("Return period",values =colsdm) +
  xlab("") + 
#!!!!! SELECT !!!!!
# ylab("Population at risk of hunger\n [% change from the no climate change level]") +
  ylab("Population at risk of hunger [% change \nfrom the median climate state]") +
#!!!!! SELECT !!!!!
  theme(
    axis.text.x = element_text(angle = 90)
  )


plot5 <- plot4 + facet_grid(RCP~Region,scales="free") 
plot(plot5)

#filenamesdmreg <- paste('../output/png/return_SDM_vsNoCC_reg.png')
filenamesdmreg <- paste('../output/png/return_SDM_vsMean_reg_N100.png')
ggsave(plot5,filename=filenamesdmreg, dpi=250,width=13, height=3.5)


#---end return period




#-----  #5 YIELD(YEXO) & price(XPRP): Boxplot across regions ----


var.in <- c("XPRP")
var.in <- c("YEXO")

for(k in 1:length(var.in)){
  varname <- var.in[k]
  
  if (varname=="PCAL"){
    labelname <- c("calorie intake [kcal/person/day]")
  } else if (varname=="PPOP_UNSH"){
    labelname <- c("population at risk of hunger [billion]")
  } else if (varname=="YEXO"){
    labelname <- c("yield [tonne/ha]")
  } else if (varname=="XPRP"){
    labelname <- c("production price [2005=1]")
  } else {
    labelname <- c("aaaaaaaa")
  }
  
#  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'_limited.png', sep = '')
  filenamebox2 <- paste('../output/png/', varname ,'/',yearname,'.png', sep = '')
  
  #colnames <- c("VAR","Year","Region","SSP","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")
  
  # Limit the num of regions
  #pyld <- pdata[pdata$VAR=="YEXO" & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD" & 
  #               pdata$Region!="XOC" & pdata$Region!="XNF" & pdata$Region!="XER" & pdata$Region!="XE25" & 
  #               pdata$Region!="USA" & pdata$Region!="TUR" & pdata$Region!="JPN" & pdata$Region!="CIS" & 
  #               pdata$Region!="CAN"& pdata$Region!="BRA" & pdata$Region!="CHN", c(-1,-2,-4)]
  
  # Limit the num of regions
  pyld <- pdata[pdata$SSP=="SSP2" & pdata$VAR==varname & pdata$Year==yearname  & pdata$SSP=="SSP2" & pdata$Region!="WLD", c(-1,-2,-4)]
  
  
  names(pyld) <- c("Region","GCM","RCP","CO2fert","SYR","N","Crop","Value","NoCC")	#Rename the column
  
  minv <- min(pyld[8],pyld[9])
  maxv <- max(pyld[8],pyld[9])
#  maxv <- c(20)
  
  #	bw <- diff(range(y$Value))/50
  #	mu <- ddply(y, "RCP", summarise, grp.mean=mean(Value))
  #	mu1 <- ddply(ypre, "RCP", summarise, grp1.mean=mean(Value))
  
  interact<-interaction(pyld$RCP, pyld$CO2fert, sep=" : ")
  
  pyld0 <- ddply(pyld,c("Region","interact","Crop"),transform,v0=min(Value))
  pyld1 <- ddply(pyld0,c("Region","interact","Crop"),transform,v25=quantile(Value,0.25))
  pyld2 <- ddply(pyld1,c("Region","interact","Crop"),transform,v50=median(Value))
  pyld3 <- ddply(pyld2,c("Region","interact","Crop"),transform,v75=quantile(Value,0.75))
  pyld4 <- ddply(pyld3,c("Region","interact","Crop"),transform,v100=max(Value))
  pyld5 <- regionlabel(pyld4)
  pyldz <- croplabel(pyld5)
  
  
  gyldbox <- ggplot() + 
    geom_boxplot(data=pyldz, aes(x=Crop, y=Value, color=interact, ymin = v0, lower = v25, middle = v50, upper = v75, ymax = v100), fill="grey80", stat="identity") +
    ylim(minv-abs(minv)*0.01, maxv+abs(maxv)*0.01) + ggtitle("") +
    ylab(labelname) + xlab("") +	MyThemeLine +
    geom_point(data=pyldz, aes(x=Crop, y=NoCC),size=6,shape=73,colour="black")
  
  gyldbox1 <- gyldbox + coord_flip()
  gyldbox2 <- gyldbox1 + facet_wrap(~Region, ncol=4)
  
  plot(gyldbox2)
  
  ggsave(gyldbox2,filename=filenamebox2, dpi=250,width=7.5, height=7.5)
  
}



#----- END --- YIELD : Boxplot across regions ----






